#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#include "gpu_util.h"
#include "cholesky_fact.h"
#include "linear_eq_band_symm_matrix.h"

#define MAX_ELEMENT 10

/*init matrix parameters based on command-line arguments and default values*/
void init_parameters(int argc, char **argv, int* matrix_order,
	int* half_bandwith, int* block_size)
{
	int arg_num = 1;
	while (arg_num < argc)
	{
		if (strcmp(argv[arg_num], "--order") == 0)
		{
			arg_num++;
			*matrix_order = atoi(argv[arg_num]);
		}
		else if (strcmp(argv[arg_num], "--hbandwith") == 0)
		{
			arg_num++;
			*half_bandwith = atoi(argv[arg_num]);
		}
		else if (strcmp(argv[arg_num], "--block") == 0)
		{
			arg_num++;
			*block_size = atoi(argv[arg_num]);
		}
		arg_num++;
	}
}

double* get_minus_one_tmp_matrix(int matrix_order, int half_bandwith)
{
	int num_elements;
	double *matrix;
	int val;
	int i, j;

	num_elements = get_band_num_elements(matrix_order, half_bandwith);
	matrix = (double*) malloc(num_elements * sizeof(double));
	for (i = 0; i < matrix_order; i++)
	{
		matrix[get_index(i, i, half_bandwith, matrix_order)] = -1;
		for (j = 1; j <= half_bandwith; j++)
		{
			val = rand() % MAX_ELEMENT;
			if (i - j >= 0)
			{
				matrix[get_index(i, i - j, half_bandwith, matrix_order)] = -1;
				matrix[get_index(i - j, i, half_bandwith, matrix_order)] = -1;

			}

		}

	}

	return matrix;
}

double* get_rand_host_matrix(int matrix_order, int half_bandwith)
{
	int num_elements;
	double *matrix;
	int val;
	int i, j;

	num_elements = get_band_num_elements(matrix_order, half_bandwith);
	matrix = (double*) malloc(num_elements * sizeof(double));
	for (i = 0; i < matrix_order; i++)
	{
		//TODO use rand for large matrices
		val = rand() % MAX_ELEMENT;
		matrix[get_index(i, i, half_bandwith, matrix_order)] = val;
		for (j = 1; j <= half_bandwith; j++)
		{
			val = rand() % MAX_ELEMENT;
			if (i - j >= 0)
			{
				matrix[get_index(i, i - j, half_bandwith, matrix_order)] = val;
				matrix[get_index(i - j, i, half_bandwith, matrix_order)] = val;

			}

		}

	}

	return matrix;
}

int calcElement(int i, int j, int matrix_order, int half_bandwith)
{
	int res = 0;

	int left = max(i - half_bandwith, j - half_bandwith);
	int right = min(i, j);

	left = max(left, 0);
	right = min(matrix_order -1, right);

	for(int k =left; k <= right; k++)
	{
		res++;

	}

	//if(i!=j)
	//{

	//	return 1;
	//}
	//return 10000;
	return res;

}

//L[i][j] = 1   returns L*L^T
double* get_positive_defined_band_matrix(int matrix_order, int half_bandwith)
{
	int num_elements;
	double *matrix;
	int val;
	int i, j;

	num_elements = get_band_num_elements(matrix_order, half_bandwith);
	matrix = (double*) malloc(num_elements * sizeof(double));
	/*if(matrix == NULL)
	{

	printf("Allocation matrix failed");
	exit(1);
	}*/

	for (i = 0; i < matrix_order; i++)
	{
		//TODO another formula for val?
		val = calcElement(i, i, matrix_order, half_bandwith);
		matrix[get_index(i, i, half_bandwith, matrix_order)] = val;
		for (j = 1; j <= half_bandwith; j++)
		{
			if (i - j >= 0)
			{
				val = calcElement(i, i - j, matrix_order, half_bandwith);
				matrix[get_index(i, i - j, half_bandwith, matrix_order)] = val;
				matrix[get_index(i - j, i, half_bandwith, matrix_order)] = val;

			}

		}

	}

	return matrix;


}

double* get_host_matrix(int matrix_order, int half_bandwith)
{
	int num_elements;
	double *matrix;
	int val;
	int i, j;

	num_elements = get_band_num_elements(matrix_order, half_bandwith);
	matrix = (double*) malloc(num_elements * sizeof(double));
	for (i = 0; i < matrix_order; i++)
	{
		//TODO use rand for large matrices
		val = rand() % MAX_ELEMENT;
		matrix[get_index(i, i, half_bandwith, matrix_order)] = (i + 1) * 100;
		for (j = 1; j <= half_bandwith; j++)
		{
			val = rand() % MAX_ELEMENT;
			if (i - j >= 0)
			{
				matrix[get_index(i, i - j, half_bandwith, matrix_order)] = (i
					+ 1) * 100 + (half_bandwith - j + 1);
				matrix[get_index(i - j, i, half_bandwith, matrix_order)] = (i
					+ 1) * 100 + (half_bandwith - j + 1);

			}

		}

	}

	return matrix;
}

void print_vector(double* v, int len)
{
	if(len > 20)
	{
		return;
	}
	printf("\nPrinting vector\n");
	for(int i = 0; i < len; i++)
	{

		printf(" %6.1f", v[i]);
	}
}

void check_solution(double * matrix, int matrix_order, int half_bandwith, double *b)
{
	for(int i = 0; i < matrix_order; i++)
	{
		int left = max(0, i - half_bandwith);
		int right =min(matrix_order -1, i + half_bandwith);
		double cur_res = 0;
		for(int j = left; j <= right; j++)
		{
			int cur_index = -1;
			if (i >= j)
			{
				cur_index= get_index(i, j, half_bandwith, matrix_order);
			}
			else
			{
				cur_index= get_index(j, i, half_bandwith, matrix_order);
			}
			
			cur_res += matrix[cur_index]*b[j];
		}
		if( fabs(cur_res - 100.) > 1e-6)
		{
				printf("\n Incorrect value of vector comp %i,  expected %i  received %f", i, 100, cur_res);
		
		}
	}
}

/*Main */
int main(int argc, char **argv)
{
	int matrix_order;
	int block_size;
	int half_bandwith;
	double* h_matrix;
	double* h_tmp_matrix; //TODO remove
	gpu_symm_band_matrix gpu_matrix;
	cublasHandle_t handle;
	long start_time, end_time, elapsed;
	double* right_hand_vector;

	init_parameters(argc, argv, &matrix_order, &half_bandwith, &block_size);
	h_matrix = get_positive_defined_band_matrix(matrix_order, half_bandwith);
	//h_matrix = get_host_matrix(matrix_order, half_bandwith);

	right_hand_vector = (double*)malloc(matrix_order*sizeof(double));
	for(int i =0 ; i < matrix_order; i++)
	{
		right_hand_vector[i] = 100;
	}

	print_band_matrix(matrix_order, half_bandwith, h_matrix);

	gpu_matrix.order = matrix_order;
	gpu_matrix.half_bandwith = half_bandwith;
	gpu_matrix.block_size = block_size;

	checkCudaErrors(cudaSetDevice(0));
	initCublasHandle(&handle);

	start_time = clock();
	copyBandMatrixToDevice(h_matrix, &gpu_matrix, handle);
	end_time = clock();
	printf("\n Time for copy %.2f seconds\n",
		(end_time - start_time + 0.) / CLOCKS_PER_SEC);

	start_time = clock();
	factorize(gpu_matrix, handle);
	end_time = clock();
	printf("\n Time for factorization %.2f seconds\n",
		(end_time - start_time + 0.) / CLOCKS_PER_SEC);

	
	copyBandMatrixFromDevice(gpu_matrix, h_matrix, handle);
	print_band_matrix(matrix_order, half_bandwith, h_matrix);

	start_time = clock();
	solve_system_factorized_matrix(h_matrix, matrix_order, half_bandwith, right_hand_vector);
	end_time = clock();
	printf("\nSystem solve time %.2f seconds\n",
		(end_time - start_time + 0.) / CLOCKS_PER_SEC);
	print_vector(right_hand_vector, matrix_order);


	h_tmp_matrix =  get_positive_defined_band_matrix(matrix_order, half_bandwith);
	check_solution(h_tmp_matrix, matrix_order, half_bandwith, right_hand_vector);

	/*h_tmp_matrix = get_minus_one_tmp_matrix(matrix_order, half_bandwith);
	copyBandMatrixFromDevice(gpu_matrix, h_tmp_matrix, handle);
	print_band_matrix(matrix_order, half_bandwith, h_tmp_matrix);*/

}
